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Abstract 

We present a method for simulating relativistic and nonrelativistic scalar field theories at finite 
density, with matter transforming in the fundamental representation of the global symmetry group 
0(N). The method avoids the problem of complex probability weights which is present in con- 
ventional path integral Monte Carlo algorithms. To verify our approach, we simulate the free and 
interacting relativistic £7(1) — 0(2) theory in 2 + 1 dimensions. We compute the two-point corre- 
lation function and charge density as a function of chemical potential in the free theory. At weak 
\4>\ 4 coupling and zero temperature we map the m? — u phase diagram and compare our numerical 
results with perturbative calculations. Finally, we compute properties of the T — fx phase diagram 
in the vicinity of the phase transition and at bare self-couplings large compared to the temperature 
and chemical potential. 
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I. INTRODUCTION 



Monte Carlo simulations provide a reliable method for studying nonperturbative physics 
of quantum field theories. It is unfortunate that with few exceptions, path integral Monte 
Carlo techniques are not suitable for the study of relativistic systems at finite density. In the 
path integral formalism, standard importance sampling techniques rely on the interpretation 
of the integrand e~ s ^ as a probability weight associated with the field configuration $. 
At nonzero chemical potential, the bosonic and fermionic Euclidean actions are typically 
complex, rendering importance sampling inapplicable. This problem is commonly known 
as the "sign problem" . 1 In many instances, the sign problem persists in the path integral 
formulation of nonrelativistic theories as well-both at finite and at zero chemical potential. 
In these cases, this sign problem does not arise because the system is at finite density, but 
rather because the action is linear in time derivatives. 

Numerous schemes have been proposed for numerical study of relativistic theories at finite 
density, with emphasis on lattice QCD. Reweighting jj] has had limited success for models 
in the regime f3fi < 1. Here, one folds the phase of the probability measure into expectation 
values of observables. Statistical ensembles are then generated using the modulus of the 
complex measure. For larger f3fi, the method fails because the calculation of expectation 
values involve intricate cancellation among phases, leading to overwhelming statistical errors. 
In addition, one must contend with the possibility that ensembles generated in simulations 
have poor overlap with the true probability distribution. In some cases, one may reduce the 
severity of the overlap problem with the use of multiparameter reweighting techniques j^]. 

Some have investigated theories at imaginary chemical potential, where lattice simulations 
are viable [l| • Using this approach, one may compute the canonical partition function for a 
system with total fixed charge Q via the Fourier transform 



In practice, this approach only postpones the sign problem since numerical integration of Eq. 
(TjQ) relies on the cancellation of phases, which become increasingly severe at large values of 
Q. A second option is to analytically continue the partition function to real \x. The analytic 

1 Some specific examples of models that avoid the sign problem include QCD at finite isospin density 
two-color QCD and low energy fermions with attractive interactions 
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continuation is typically based upon a Taylor series expansion about /x = and therefore 
this approach is only valid within the radius of convergence-inevitably failing at and beyond 
any critical point. In addition, the analytic continuation is limited by the periodicity of the 
partition function Z(ifi) = Z(ifi + i2k/ f3) along the imaginary /j axis-a property which is 
evident from analysis of the fugacity expansion. 

In the case of scalar theories, which is the focus of this paper, we investigate the prospects 
of finding alternative representations for the partition function which are suitable for Monte 
Carlo simulations. To this end, we formulate relativistic and nonrelativistic O(N) scalar 
theories io terms of dual lattice variables Q. The method makes use of U(l) character 
expansions to explicitly integrate out the angular mode of a complex scalar field. In the 
dual representation, the angular mode is replaced by a set of integer valued link variables 
which may be interpreted as the conserved current density of the system. The path inte- 
gral formulation is particularly suited for Monte Carlo simulation because complex phases 
associated with the chemical potential (in the case of relativistic theories) and lattice time 
derivative (in the case of nonrelativistic theories) are avoided. All physical observables are 
measurable within the formalism-free of sign problems-provided the corresponding opera- 
tors are neutral in charge. We will consider scalar theories in the grand canonical ensemble, 
although as we shall see, the formulation is applicable to canonical ensembles as well. 

To motivate the approach, we begin by considering a quantum mechanical rigid rotor at 
finite temperature (T = (3~ l ) and chemical potential /i coupled to angular momentum. This 
model corresponds to a (0 + 1) -dimensional relativistic scalar theory at finite charge density, 
where the dynamics of the heavy radial mode have been ignored (or integrated out). The 
continuum partition function is defined by the path integral representation 




where S(r) = e %e ^ is an element of U(l) and I is the moment of inertia. At finite tem- 
perature, we impose periodic boundary conditions on £ in the imaginary-time direction. 
The periodicity of E in the imaginary-time direction translates into the boundary condition 
9(r) = 9{t + (3) + 2ttw on the angular mode, where w is the winding number associated with 
the field configuration 9{r). Our definition of the partition function implicitly includes a 
sum over all values of the winding number; this sum may be made explicit by introducing a 
new field 9 w (t), which represents configurations with fixed winding number to. The partition 
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function may then be written as: 



(3) 



By taking advantage of the correspondence between configurations with zero and nonzero 
winding number, namely 9 w (r) = Oq(t) + 2ttwt/i3, we may express the partition function as 
the product of two terms: 



where the chemical potential decouples from 6q but remains coupled to winding number. 
Since we are only interested in the /i-dependence of the partition function, we may neglect 
the /3-dependent path integral over the zero winding number sector of the theory. Using 
the Poisson summation formula, the partition function is expressed as a sum over positive 
weights: 



The resulting representation for the partition function is given in terms of the momentum 
conjugate variable to the coordinate 9, namely the angular momentum Q, indicating that the 
sign problem is a basis dependent problem; this observation forms the basis for formulating 
the path integral for scalar field theories in terms of dual lattice variables. In passing from 
Eq. (jlj) to Eq. (jSj), the dimensionless ratio T/(3 is mapped to (3/X-a, typical property of 
duality transformations-and the chemical potential now couples to total angular momentum 
(or charge) Q of the system. 

In Sec. ([IT]) we will develop a generalized formalism for relativistic and nonrelativistic 
O(N) scalar field theories defined on a space-time lattice and describe how to compute 
expectation values within this framework. In addition, we outline standard procedures for 
solving certain continuity constraints that arise in our derivation. For completeness, we 
show how U(l) gauge fields fit into the formalism as well. In Sec. fillip we explain several 
numerical algorithms for simulating relativistic bosonic theories at finite density and discuss 
additional issues that require resolution in order to render such simulations practical. Sec. 
(HVj) pertains to the details of our simulations. A comparison of numerical and analytic 
results is presented for the free and interacting |0| 4 theory in 2 + 1 dimensions and at weak 
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coupling. Sec. (jVj is devoted to an application of our methods to the same model at finite 
temperature and density, and in a nonperturbative regime. Here, we study the behavior of 
the theory in the vicinity of the zero temperature and density critical line, and determine 
several universal constants associated with the phase transition. 

II. FORMALISM 

For simplicity, we consider a single complex field <fi with chemical potential fi coupled 
to the U(l) current associated with the global transformation — > e ta <j). Generalizing the 
following results to the case of O(N) with fundamental matter is straight-forward, since 
the chemical potential always couples to U(l) charges embedded within O(N). If N is 
even, then with an appropriate change of basis, the theory can be written in terms of N/2 
complex scalar fields — each of which is coupled to an independent chemical potential. The 
techniques discussed for a single complex scalar field can then be applied to each of the N/2 
fields individually. If iV is odd, then there will be an additional real, uncharged scalar field 
which plays no role in the following discussion. 

For a single complex scalar field, the grand canonical partition function in Euclidean 
space-time takes the generic form 



where So and Si correspond to the free and interacting parts of the action respectively. The 
specific form of So will depend on whether one is interested in relativistic or nonrelativistic 
scalar theories; both cases will be discussed in detail in the following subsections. We focus 
primarily on the free part of the action at nonzero /z, though interactions and additional 
flavors may be added in a trivial manner. We assume that Si is symmetric under U(l) gauge 
transformations for reasons which will become clear later on. 

We regulate the (d + l)-dimensional continuum theory by applying the usual lattice dis- 
cretization procedure. We set the spatial and temporal lattice constants equal and work on 
a Nf x N T lattice with an inverse temperature given by (3 = N T . For notational convenience, 
all dimensionful parameters are written in units of the lattice constant. Continuum deriva- 
tives are replaced with finite differences d u <p — > <Pn — 0n-e„, and following the prescription 
of 0], we introduce the chemical potential as the imaginary zero component of the vector 
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potential (i.e. A v — > i{j,8 U} o). This prescription not only avoids spurious cutoff dependence 
in thermodynamic observables, but is also a natural choice within our formalism, enabling 
us to draw a direct connection with the fugacity expansion 

Throughout the paper we impose periodic boundary conditions in the space and time direc- 
tions. 

We begin by applying the dual lattice transformation to the partition function for free 
relativistic and nonrelativistic scalars. Later, we introduce interactions as well as repeat 
the derivation for expectation values of observables. A general discussion on character 
expansions can be found in 0, Q] . 



A. Relativistic theory 



The free relativistic Euclidean lattice action is defined by 



n Vn l 
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where v runs over positively oriented basis vectors. In terms of the polar coordinates 
pe* e , the partition function may be expressed as 
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As previously stated, performing the ^-integration analytically will yield a new represen- 
tations for the partition function, expressed as a sum over integer valued link variables as 
well as the radial degrees of freedom. To facilitate the integration over angular degrees of 
freedom, we make use of the convergent series expansion 
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where Ii(x) is the modified Bessel function of order I, with i6R and z G C. Inserting Eq. 
f fTUj) into the partition function and exchanging the sum over modified Bessel functions with 
the product over nearest neighbors, we obtain 

-2tt 



ZM= [d6] / [pdp] £ J] e-( 2d+2+ - 2 )^n^.( 2 ^^+eJe«- 
{in,v}ez n L v 
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where the sum over integer valued link variables I correspond to the sum over characters of 
U(l) (or Fourier modes). For the moment, the sum over link variables in equation Eq. ffTT]) 
is unconstrained. However, upon integration over 9 one finds that the only nonvanishing 
contributions to the partition function are those terms independent of 9. These contributions 
may be characterized by link field configurations that satisfy the continuity equation 

d ■ l n = J2 + = > ( 12 ) 

V 

where / n i/ = —l n +e v -u- The partition function takes the final form 

Zr(v)= T / W"' 5 *^^. (13) 



d-l n =0 



where 
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(14) 



The constraint given by Eq. ffl2|) allows us to identify the link field I as a conserved current 
associated with the U(l) symmetry of this theory. This observation will be justified later 
when we introduce gauge interactions. Furthermore, the chemical potential now couples to 
the the total charge Q = J2 n ln,o- 

A key property of the modified Bessel functions appearing in Eq. (fT4"|) is that Ii(x) > for 
all x > and / e Z. The partition function is therefore a sum over real and positive weights, 
rendering the theory suitable for Monte Carlo simulations. We would like to emphasize at 
this stage that our derivation remains unchanged in the presence of U(l) gauge symmetric 
interactions because such interactions are necessarily independent of the angular mode 9. 



B. Nonrelativistic theory 



In the nonrelativistic limit, one may proceed in an analogous fashion to the relativistic 
case. 2 The free Euclidean lattice action is defined by 
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(15) 



2 A similar approach was taken by [9( for the Bose Hubbard model, starting from the operator formalism. 
Numerical simulations were performed in the so-called "hard core" approximation. 
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where the sum over basis vectors run over the spatial directions. After performing shifts 
in the dummy variable n and changing to polar coordinates, the action takes the form 
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Using Eq. (TTDIl and the relation 
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we integrate out the angular degrees of freedom and obtain the dual partition function 
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with the nonrelativistic dual action defined by 

ol (Pnpn+e ) ln ° 



Snr[p,1] 



n 



-(l+d/m)p2 _ 



(19) 



As with the relativistic theory, the nonrelativistic current density is constrained by the 
continuity equation Eq. (|T2l) . In addition, the zero-component of the current density is 
restricted to non-negative integers, reflecting the fact that nonrelativistic particles may only 
propagate forward in time. Once more, the partition function is free of sign problems and 
is therefore suitable for Monte Carlo studies. 



C. Observables 

We next address the computation of expectation values for both the relativistic and 
nonrelativistic theories. We begin by considering the generic 2fc-point correlation function 

G mtLm fc ,n lv ..n fe = (C • • • tmrfm • • • 

= (Pm a • • .Pm fc p ni • .. Pak e i ^-^ k -e mi -...-e mk ) ) (20) 

Following the formalism in the previous subsections, we integrate over angular variables and 
obtain 

1 f°° 

Ggt,™ k ,ni,-..n k =z7-S E / [prfp]p mi ---Pm fe Pn 1 ...p nfc e- S M e ME„^C ) (2l) 

d-l n = J n 
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where S may represent either the relativistic or nonrelativistic action, including the inter- 
action term Si. The sum over current densities is, however, no longer divergence free. In 
heuristic terms, an insertion of <f) n creates one unit of charge at site n, whereas an insertion 
of <p* m annihilates a unit of charge at site m. We therefore expect the nontrivial divergence 
condition 

■ l n — <5n,ni + • • • + $n,n k ^n,mi • • • ^n,m fe = Jn > (^2) 

where J represents a linear combination of sources and sinks. In general terms, all gauge 
variant operators will give rise to nontrivial divergence properties of the current density. We 
discuss solutions to these constraints in the following section. 

Applying the formalism above, it is possible to calculate all thermodynamic observables of 
interest. The dual representation is particularly well-suited for calculating current densities, 
charge distributions and current density correlation functions, since the current density is 
one of the dynamical degrees of freedom in this representation. In our numerical studies of 
relativistic bosons, we pay particular attention to the charge density (n) and susceptibilities 
Xab, defined by 

Xab — N T N% ((AB) — (A)(B)) , (23) 

where a and b represent source terms for the generic operators A and B. For this study, 
we consider the susceptibilities X^m 2 an d Xm 2 m 2 where the source terms fi and m 2 
correspond to the operators 

A n 

and 



T^s n 



m N T N? 

respectively. 



^E/^ ( 25 ) 



D. Constraints 

Standard procedures exists for solving the divergence constraints on current densities {f| . 
We begin by decomposing the link variables in terms of homogeneous and particular solutions 
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FIG. 1: A particular solution to d ■ l n = J n arising from the four-point correlation function 
Bold link variables equal +1( — 1) if the link orientation and basis vector are parallel (antiparallel) . 
All remaining links equal zero. 

to d ■ l n = J n . The homogeneous solutions may me labeled by 1® and are characterized by 
d + 1 integer valued topological charges defined by 

= • (26) 

n 

As previously discussed, Qq corresponds to the total conserved charge of the system, whereas 
Qj corresponds to the current flux in spatial directions. With this decomposition, expecta- 
tion values take the generic form 

d-l%=0 

Note that the particular solution l J is not unique; one solution may be related to another 
by a suitable change of variables. A convenient choice for l J is to set all links equal to zero 
except for a set of "strings" which connect source-sink pairs as shown in Fig. (CQ). In the 
relativistic theory, such strings may flow either forward or backward in time, and may wrap 
completely around the space-time torus; this would correspond to an increase or decrease 
in the total charge of the system by one unit (note that this may be undesirable in the 
canonical ensemble). In the case of nonrelativistic theories, we shall assume for simplicity 
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FIG. 2: Solving zero divergence constraints by introducing integer valued plaquette variables. The 
sum over all (oriented) plaquettes P circulating a common link defines the value of link variable I. 

that such strings only flow forward in time and do not wrap entirely around the space-time 
torus. 

The zero divergence constraint on on link field configurations in the Q = sector may be 
solved by changing to plaquette variables P n>fiu = —P n ,u^ G Z. A plaquette is specified by a 
site n and the two (different) basis vectors /i and v that span it and is positively oriented if 
< v. The divergence free link variables in the Q = sector are given by 



which is depicted graphically in Fig. ([2]). Note that Eq. (1281) has a continuum analog; in 
three dimensions, for instance, a divergence-free vector field may always be written as the 
curl of an unconstrained vector field. Similar relations hold in higher dimensions as well, in 
the language of differential forms. Configurations associated with the Q ^ sectors may be 
obtained from the Q = configuration by increasing the value of all link fields belonging to 
a closed loop around the space-time torus by Q u units in the v direction. 

E. U(l) gauge fields 

In the presence of gauge interactions, the sign problem reemerges within our formalism. 
In the compact U(l) gauge theory (scalar electrodynamics), the vector potential A couples 




(28) 
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to the integer valued conserved current I. This result is easily verified by repeating the 
derivation in the preceding subsections, which yields the gauge invariant interaction 

SDJ^^InA- (29) 

n v 

Note that the form of the gauge coupling confirms our interpretation of the link field as a 
current density. As one may expect, by going to the dual formulation we have traded the 
chemical potential sign problem for a vector potential sign problem. A natural remedy is 
to apply the character expansion to the gauge partition function in addition to the scalar 
partition function. Such transformations are not new in the context of gauge theories; 
character expansions have played an important role in understanding the strong coupling 
behavior of such theories. In the case of lattice Yang-Mills theories with Wilson-type actions, 

n 

the character expansion has a finite radius of convergence [10|], which likely limits the utility 
of this discussion. 

For completeness, we briefly describe the application of character expansions to dynamical 
U(l) gauge fields coupled to the conserved current I. Consider the compact U(l) gauge 
theory defined by the Wilson action 

Z gauge = \ [dA] 6^ En E ^ ^ En E " , (30) 

JO 

where Q n<f iu = A nfl + A n+eM>i/ — A n+ei/fl — A n v is the lattice field strength tensor. The gauge 
partition function vanishes unless the topological charge Q v associated with the conserved 
current I equals zero. After applying the character expansion to the gauge partition function, 
we obtain 

Z gauge = £ ]J ]J I Fnjuf (^) , (31) 

where the integer valued plaquette variables F satisfy 

{p ' F-a) p = ^ {F n ^ u + i^ n ,/i— u) = ln,u • (32) 
v 

This constraint is simply the lattice analog of the continuum equation of motion for the 
gauge field. Note that the partition function Eq. (13T]) implicitly depends on the current 
density / though the Eq. 
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In the absence of current sources, there are d(d— l)/2 topological charges associated with 
the plaquette variables defined by 

(33) 

n 

These charges may be visualized as two-dimensional oriented "sheets" that wrap around the 
space-time torus. The source-free constraint on the plaquette field F may then be solved by 
a change of variables to integer valued cubes. Inhomogeneous solutions to d ■ F n = l n may 
be obtained by joining all current loop source-sink pairs with flux tubes, in an analogous 
manner to our treatment of current sources and sinks. 

III. ALGORITHMS 
A. Loop algorithm 

The dual partition function for scalar theories involves an integral over site variables 
p and a constrained sum over the link field I. As discussed in the previous section, the 
constraint on I is solved by making an appropriate change of variables. For the purpose 
of numerical simulation one need not work directly with plaquette variables, however. An 
easier and equivalent approach is to simply consider fluctuations in conserved current loops. 
We consider two types of current loop updates: local loops shown in the lower left of Fig. 
([3]), and global current loops shown in the upper right of Fig. ([3]) which correspond to an 
increase or decrease in the charge of the system by one unit. Note that an algorithm based 
on local and global loop updating is ergodic: any divergence-free configuration with charge 
Q may be obtained from any other with like charge by a sequence of local current loop 
updates. Divergence-free configurations in different topological sectors may be related by 
global current loop updates, possibly succeeded by a sequence of local loop updates. 

One may use the standard Metropolis [11] accept/reject method for local site and loop, 
as well as global loop updates. However, as lattice volumes become large, the acceptance 
rate for the global loop update is expected to drop to unacceptably low levels. This in 
turn adversely affects the correlation time for configurations within the ensemble. The low 
acceptance rate is due to the fact that the charge transition probability scales as e~ L , where 
L ~ N T , N s is the linear length of the global current loop. The likelihood for the system 
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FIG. 3: Local (lower left) and global (upper right) current loops on a periodic lattice. Bold links 
are updated by +1(-1) if the link orientation and basis vector are parallel (antiparallel). 

to become stuck in charge sectors which are energetically disfavored increases with volume 
sizes; this issue may be particularly problematic during the equilibration stage of a Monte 
Carlo simulation. 

The problem of slow tunneling rates may be handled in several different ways. In the 
the thermodynamic limit, fluctuations in the charges Qj may be neglected altogether since 
these fluctuations will have little influence on local observables. Qo fluctuations may be 
neglected for similar reasons at zero temperature, as well as in cases where canonical ensem- 
ble simulations are of interest. At finite but large temperatures and in the grand canonical 
ensemble, one might avoid the problem by implementing an intermediate equilibration step 
to facilitate fluctuations between different Qo sectors, followed by an accept/reject step. 



B. Worm algorithm 

An alternate method for avoiding long correlation times associated with charge fluctua- 
tions is the "worm" algorithm [12 ] . This algorithm may be directly applied to the partition 
functions given by Eq. (|T3l) and Eq. (118p which were derived from the character expansion 
in addition to partition functions derived from hopping parameter expansions. We will focus 
on the former; the later is discussed in [12] at zero chemical potential and may be generalized 
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to dense matter in a straight-forward manner. We emphasis that in each of the above cases 
the role played by the chemical potential is the same. 

The worm algorithm employs a dynamical source and sink associated with the current 
density to generate link field configurations. The process begins by associating a current 
source and sink with a randomly chosen site on the lattice. One then attempts to move 
the source to a neighboring site by randomly updating one of the 2{d + 1) link fields (us- 
ing Metropolis accept/reject method) associated with the bond by one unit. As links are 
updated, the divergence condition on the link field becomes d ■ l n = 5 n na — 5 njn6 , where n a 
and rib represent the respective locations of the source and sink. Proceeding one link at a 
time, the source traces out a string or worm on the space-time lattice. Most of the link 
field configurations generated by this method are unphysical and may be discarded because 
they do not satisfy the divergence- free constraint. If the endpoints of the string occupy the 
same site, then the divergence-free constraint is satisfied and the field configuration is stored 
for later calculation of physical observables. A new site is then chosen at random and the 
process is repeated. 

At finite density, the probability for the source to propagate forward in imaginary time 
is enhanced by a factor of e M . This bias magnifies the likelihood for the source to traverse 
the temporal extent of the lattice after a sequence of updates, resulting in an increase or 
decrease of the topological charge Q by one unit. The main advantage of this algorithm is 
that tunneling processes are avoided. Instead of going over the energy barrier associated with 
changes in Q , the worm algorithm allows one to circumnavigate the barrier by accessing 
unphysical states of the system-states where the charge Q is not well-defined. 

C. The gauge sector 

In numerical simulations involving both gauge and matter fields, one must allow fluctu- 
ations in link variables / and plaquette variables F, subject to constraints. The issues may 
be dealt with in an analogous fashion to the loop algorithm. The divergence constraint on 
the plaquette field is satisfied by considering fluctuations in local flux tubes (cubes). On 
small lattices, fluctuations in the global charge Q^ v must also be considered. Since these 
charges correspond to sheets (with area A), the probability for fluctuations of this kind 
scale as e~ A . In addition to these, one must consider fluctuations involving gauge plaque- 
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FIG. 4: Local flux tube (left) and flux tube connecting a current loop source and sink (upper 
right). 

ttes bounded by current loops; these current loops act as flux tube sources and sinks. An 
example configuration is displayed in Fig. (j3J). 



D. Reweighting 

Before efficient simulations can be performed, additional issues must be addressed. Ideally 
one would like to generate a single set of configurations {$j} = {p i: where i — 1, . . . , Af, 
and M is the sample size, distributed according to the probability weight e -5 '*'. From this 
ensemble, one would then like to approximate the expectation values of all observables using 
the relation 

i=i ^ ' 

Unfortunately, in the dual representation link variables must satisfy specific divergence con- 
straints which are governed by the gauge transformation properties of each observable. As 
a result, configurations generated in a simulation may only be used to calculate expectation 
values of those observables which transform similarly under gauge transformations. With 



the use of reweighting J4], this issue may be avoided. Let (0(p))j represent the expectation 
value of some p-dependent operator in the presence of source J. Then referring to Eq. (J2j 
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one can show that 



(0( P ))j 



0(p)e~ §ip ' lQ+lJ)+ ^ l i° 



(35) 



e -S( P ,lQ) 



o 



where (. . .)o is the expectation values with respect to source- free configurations. In this 
scheme, all expectation values may be computed from the same source-free distribution. 
One disadvantage of reweighting, however, is that one must now evaluate expectation values 
of nonlocal observables. Reweighting is expected to give reliable results provided the Monte 
Carlo distribution and target distribution have sufficient overlap. In particular, the extent 
of correlation functions should remain well within the bulk of the lattice. Otherwise, the 
Monte Carlo and true distributions may lie over different topological sectors, leading to large 
systematic errors. One may be able to reduce the systematic errors arising from the overlap 
problem by taking advantage of the fact the the particular solution l J is not unique. In 
particular, one could average the left hand side of Eq. fl35|) over several different particular 
solutions l J , thereby improving estimates of the expectation value. While numerical results 
will most likely improve, the systematic errors remain uncontrolled with this scheme. 

IV. SIMULATIONS AT ZERO AND WEAK COUPLING 

We have now provided a framework for studying scalar theories at finite chemical potential 
on the lattice, free of sign problems. We test the convergence of our algorithms by performing 
Monte Carlo simulations of the (2 + l)-dimensional free, relativistic U(l) theory at finite 
density and finite temperature. The algorithm is then applied to the same theory with an 
additional repulsive interaction 



Here, we determine the m? — fi phase diagram at infinite volume and zero temperature, and 
in a regime where perturbation theory is reliable. 

Since these simulations are exploratory in nature, we work in 2 + 1 dimensions and limit 
our lattice sizes to 10 2 x 10 or smaller, employing the loop algorithm discussed above. During 
each run, we update the system through ~ 2 x 10 6 Monte Carlo steps (1 Monte Carlo step 
= 1 sweep through the system) and estimate that the correlation times are on the order 
of r corr ~ 2 x 10 2 steps. Following an equilibration time r eq ~ 2 x 10 4 , configurations are 




(36) 
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accepted every r corT steps to ensure that ensembles are uncorrelated. We include fluctuations 
in global charges Q v to ease the comparison of numerical results with analytic calculations. 



A. Free theory 

We begin by measuring properties of the free theory at finite chemical potential. The 
free theory is exactly soluble because the action is a quadratic form involving <fi. At zero 
coupling and finite lattice spacing, the theory is stable provided 4sinh 2 (/i/2) < m 2 . At finite 
volume, the free energy of the system is given by 

-N T N s f {n) _ JJ 



e 
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77^ + 4 sin ' 



p 1 V 2 / 

11 r^vh ( AT pn C |,-lf A _ nnah I AT ,A ' 



cosh(A T cosh 1 £ p ) — cosh(iV T /x) 

where 




£ P = l + ± U 2 + 4£sin 2 (f) , (38) 



and the momenta take the discrete values p u = 2im u /N u (Nq = N T and Nj = N s ) within 
the Brillouin zone. From the free energy we compute the charge density as a function of /a, 



(n) 



i-y shlh } N ^ . (39) 

N* ^ cosh(A r cosh" 1 ^ p ) - cosh(A r /i) 



In the low temperature regime (/3m ^> 1) the charge density remains zero for 4sinh 2 (/i/2) < 
m 2 . At higher temperatures (/3m ~ 1) particles (antiparticles) become thermally excited, 
leading to finite densities for |/i| < m. In both 4sinh 2 (/i/2) — > m 2 the charge 

density n diverges, signaling the onset of Bose-Einstein condensation. 

Simulations of the free theory are performed at a fixed mass (m = 1) and for chemical 
potentials ranging from li = 0.0 — 0.9. Fig. (jSJ) provides a comparison of the exact analytic 
and numerical computation of the charge density as a function of chemical potential for 
(3m = 10. In addition, we measure the two-point correlation function = (0n+re o ^n) 
using the reweighting procedure described in Sec. fillip . Our numerical results for in 
Fig. (J5]) agree with exact calculations for time separations r <C N T . For time separations 
t > N T , the reweighting method fails to yield correct results. This is presumably due to 
poor overlap between Monte Carlo and target distributions. 

18 



0.02 
0.015 
<n> 0.01 
0.005 


0.5 0.6 0.7 0.8 0.9 

FIG. 5: Charge density (n) as a function of fi for the free relativistic theory (m = 1) on an 10 2 x 10 
lattice. Dashed line indicates the exact analytic calculation, data points indicate numerical results. 

B. Interacting theory 

At weak coupling (A = 1) and zero temperature we determine the critical mass m 2 
as a function of //. The transition occurs when the effective mass m 2 jy(m 2 ,/i) = m? — 
4sinh 2 (/i/2) + S(m 2 ,/i) vanishes, where S(m 2 ,/i) represents the selfenergy correction to the 
free propagator. We use this relation to compute m 2 for fixed values of fi and A. At tree 
level, m 2 is simply given by m 2 = 4 sinh 2 (/z/2); this line is indicated by the dashed curve in 
Fig- ©• We compute the selfenergy perturbatively in the infinite volume limit and at finite 
lattice spacing. The first-order correction arising from the 1-loop diagram shown in Fig. (JHJ) 
yields the relation m 2 = 4sinh 2 (/i/2) — T^i^m 2 ,^). The 1-loop self energy is proportional to 
the free energy f Q of the noninteracting theory and in the infinite volume limit one finds 

Ei(, "'" )=A L(lvrT i '( 1 - eM (^v^ 1 )) • (40) 

where £ p is given by Eq. fl38|) and H(x) represents the Heaviside step function. The domain 
of integration is over spatial momenta pj G [— 7r,7r]. We solve for m 2 to leading order in A 
by substituting our tree level result for m 2 into E^m 2 ,/!). The first order correction to the 
critical line is indicated by the solid curve in Fig. ([7]). Above the critical line (m > m c ) 
the charge density is strictly zero at zero temperature. Below the critical line (m < m c ) the 
charge density is positive for fi > and negative for /i < 0. 

In our numerical simulations, the critical mass is determined by observing the dependence 
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FIG. 6: Two point correlation function G\ for the free relativistic scalar theory with /i/m = 0.0 
and 0.9. Error bars indicate only the statistical uncertainty and not that due to systematic errors. 
Exact analytic calculations (squares) agree with Monte Carlo simulations (triangles) for r < N T /2. 
Results were obtained for a 10 2 x 10 lattice Discrepancies between the numerical data and analytic 
results for r > N T /2 are discussed in the text. 

of the susceptibilities Xw an d X^m? on the linear size N T , N s of our lattice. In the infinite 
volume limit, these quantities are discontinuous (assuming \x ^ 0) across the phase boundary. 
At finite volume, however, it is well known that no phase transition exists; discontinuities in 
thermodynamic quantities become smooth cross-overs. We obtain an approximation to the 
infinite volume critical mass by plotting Xnn anc l X^m 2 as functions of m 2 for various lattice 
sizes. A typical set of data is displayed in Fig. Q for \i = 0.4. Simulations are performed 
at a trial mass m 2 , near m 2 c . Susceptibilities in the neighborhood of m 2 , are then obtained 
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FIG. 7: The charge density phase diagram in the (//,m 2 ) plane at A = 1. The dashed (solid) 
curve corresponds to the tree level (1-loop) calculation of m 2 as a function of fj,. The upper (lower) 
region corresponds to the unbroken (broken) phase. Data points are obtained from the volume 
dependence of Xum 2 i as shown in Fig. ([9]). 



In the limit N T , N s — > oo, the susceptibility curves should intersect at the same point, 
revealing the infinite volume value of the critical mass. In practice, one should apply a 
finite scaling analysis to accurately determine m 2 c in the thermodynamic limit. Since our 
lattices are relatively small (scaling analysis is only valid for N T , N s ^> 1) and our computing 
resources limited, we choose not to do so here. Estimates of m 2 c are obtained from the 
N T , N s = 4, 6, 8, and 10 crossings. In the special case /i = 0, the average number density and 
susceptibilities remains zero across the phase boundary. Here, we study the the properties 
of the charge distribution associated with Q as a function of lattice volume. In particular, 




FIG. 8: ldoop self energy diagram Si(m 2 ,/x). 



from the same ensemble with the use of 



(0) m 2 




(41) 
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FIG. 9: Volume dependence of x^m 2 f° r A 4 = 0.4 and lattice sizes ranging from N T ,N S = 4 — 10 
in the space and time direction. Curves are obtained by extrapolating data from m 2 = —0.057 
(dashed line) via reweighting. 

in the infinite volume limit the Binder cumulant U = 1 — (Qq) /3(Qq) 2 equals zero in the 
broken phase (m 2 < m 2 ) and diverges in the unbroken phase (m 2 > m 2 ). Our numerical 
results for the m 2 — /i phase diagram shown in Fig. (JTj) appear consistent with perturbation 
theory at weak coupling. 



V. AN APPLICATION 



We now turn to an application, focusing for simplicity on the relativistic U(l) lattice 
theory with the interaction given by Eq. fl36|) and in a nonperturbative regime. The pur- 
pose of this section is to show that full scale simulations are feasible using the methods 
outlined above. Since our path integral representation is derived in the imaginary-time for- 
malism, the utility of our approach is limited to systems in thermodynamic equilibrium. 
Nonetheless, the algorithms may, for instance, be used to study the universal properties of 
the (2 + l)-dimensional relativistic |0| 4 theory, which is known to exhibit a second order 
phase transition at zero temperature and density. At (3, V = oo and at the phase transition, 
the theory exhibits no intrinsic scale. As a result, one can make predictions about the de- 
pendence of observables on the physical parameters of the theory based on universality and 



the scaling hypothesis. As an example, dimensional analysis suggests that (n) = Bfi 2 [1 
and T c = C\i I4J] as one moves away from the zero temperature and density critical line 
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FIG. 10: to 2 dependence of U for A = 192, fi = 0.0 and space-time volumes N T , N s = 20 — 36. 

The proportionality constants B and C are believed to be universal and nonperturbative, 
and therefore ideal for calculation by numerical means. 

For the calculation of universal constants, we choose the self coupling A = 192 and 
determine the critical mass by studying the volume dependence of the Binder cumulant 
U(m 2 ) at zero chemical potential. Results for U(m 2 ) were obtained using the worm algorithm 
in order to facilitate charge fluctuations at larger volume sizes, as discussed in Sec. (llllj) . 
Scans in the mass parameter were performed for lattice sizes ranging from N T , N s = 20 — 36 
sites. The results plotted in Fig. ( fTOl) suggest that the critical mass is approximately 
m 2 ~ —26.05 for this choice of A. 

The universal constant B is determined by studying the volume dependence of the ratio 
(n)/fi 2 at A = 192 and m 2 = m 2 . Curves were obtained by direct measurement of the 
number density, once again using the worm algorithm. Scans displayed in Fig. ffTTj) were 
performed for values of the chemical potential ranging from \i = 0.1 — 0.5 and lattice sizes 
ranging from N T ,N S = 20 — 36 sites. In the infinite volume limit, (n)/fi 2 is expected to 



converge to B. Previous indirect measurements of this quantity have found B = 32(1) [15 ]. 
With exception to the \x = 0.1 level curve, we find B ~ 0.32 for large space-time volumes-in 
apparent agreement with the earlier calculation. 

Turning now to the finite temperature behavior near the critical line, we study the tem- 
perature dependence of the number density for chemical potentials ranging from /i = 0.1—0.5 
and using the loop algorithm described in Sec. fllllj) . Curves displayed in Fig. (112j) were ob- 
tained by direct measurement of the number density at a spatial volume N s = 16. Near the 
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FIG. 11: Volume dependence of (n) / for A = 192, m 2 = —26.05 and yt, = 0.1 — 0.5. 

critical point (A,m^) = (192,-26.05), the effective theory describing the Goldstone mode 
9 predicts weak temperature dependence in the T <C fi regime. More precisely, using the 
arguments of the leading order contribution to the continuum effective Lagrangian is 



B 



213/2 



£= 3 [(d T 6 + + (V6) 2 ] 
Expanding the Lagrangian in inverse powers of //, we find 



{d T 9f + -{V9f 



+ 



and to leading order in temperature, the number density is given by: 



(42) 



(43) 



(44) 



This predicted temperature dependence is evident in our numerical calculation of (n) j \i as 
a function of T for T fi. Results of this calculation are shown in Fig. fjl2H . 

For a precise determination of the universal constant C, which characterizes the shift in 
T c as a function of fi, we fix the inverse temperature at iY> = 4, 6 and 8 and study the 
//-dependence of the susceptibility Xm 2 m 2 f° r spatial lattice volumes N s = 24, 36 and 48. 
Results are shown in Fig. f|T3|) . In the thermodynamic limit iV s — > oo, this susceptibility is 
expected to diverge as one approaches criticality. From the infinite volume location of the 
maximum /i maa: (A^ s — > oo) of the susceptibility, one can determine the universal constant 
by using the relation C~ l = /3fi max (N s — ► oo). We leave an infinite volume extrapolation of 
Hmax{N s — > oo) to a later study; at finite volume we find the approximate value of C ~ 0.44 
for each choice of (3 and for the spatial volume N s = 48. 
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FIG. 13: jJL dependence of Xm 2 m 2 f° r ^ = 192, and m 2 = —26.05 for spatial volumes N s = 24 
(diamond), N s = 36 (star) and N s = 48 (square). 

VI. CONCLUSION 

We have derived a new representation for the bosonic partition function which may be 
interpreted as a path integral over current densities. The representation avoids the problem 
of complex phases and is therefore suitable for use in numerical simulations. The formalism 
is applicable to relativistic and nonrelativistic O(N) scalar theories at finite density and in 
the case N = 2 allows for the inclusion of gauge interactions at strong coupling. We have 
verified the method by performing numerical simulations of the free relativistic scalar field 
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theory in 2 + 1 dimensions, finding that numerical results for the two-point correlator and 
charge density agree with exact calculations. The method was then applied to |0| 4 theory at 
weak coupling where we have correctly reproduced the m 2 — /i phase diagram. Finally, we 
have demonstrated the utility of the formulation by performing numerical calculations of the 
universal constants B and C associated with the phase transition at finite temperature and 
chemical potential. We find the approximate values B ~ 0.32 and C ~ 0.44; a full analysis, 
including infinite volume extrapolation of the data is left for future studies. 

The methods and ideas discussed in this paper may be extended to include the addition 
of higher representation fields and multiple flavors. In cases such as this, positivity of the 
dual partition function will depend on the specific form of interactions and whether of not 
the resulting character expansion coefficients satisfy the desired positivity requirements. 

Finally, we have been unable to generalize our method to the case of fermionic systems 
for the obvious reason: no polar decomposition exists for Grassmann variables. Similar 
methods which are based on the hopping parameter expansion are known to fail as well 
because configurations with an an odd number of fermion loops in the expansion carry 
negative weight. It is likely that the fermion sign problem is fundamentally different from 
that for bosons and the techniques describes here are simply inapplicable. 
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